#include "Podzial.h"
#include <cmath>
#include "FilesManager.h"
#include <iostream>
#include <iomanip>

#ifdef GPU_INCLUDED
        #include "GPGPU.h"
#endif


/* ----> METODA XY - DLA KSZTALTOW ZBLIZONYCH DO PROSTOKATA. */
/* ----> procedura wypelnia wnetrze sieci. */
/* ----> xnew(j) - wspolrzedne j-tego wezla tworzonej sieci: */
/* ---->         xnew(2*j-1) - wspolrzedna x-owa */

/* ---->         xnew(2*j)   - wspolrzedna y-owa */


//GLOB.nkp.mw, GLOB.nkp.nw, GLOB.zargd.xpodz,
//GLOB.zargd.ypodz, &GLOB.nkp.nx, GLOB.zargd.xp, GLOB.zargd.yp,
//GLOB.zargd.rp, GLOB.zargd.xrp, GLOB.zargd.yrp, GLOB.zargd.ttp,
//GLOB.zargd.alpp, GLOB.zargd.alkp, GLOB.zargd.xw,
//&GLOB.nkp.xcharl, &GLOB.nkp.ycharl, &GLOB.nkp.xchar, &GLOB.nkp.ychar,
//&GLOB.nkp.xchard, &GLOB.nkp.ychard, &GLOB.nkp.xcharld,
//&GLOB.nkp.ycharld, typ

//#define GPU_INCLUDED

int Podzial::podzial(int mw, int nw, double *xpodz,
                     double *ypodz, int *nn, double *xn, double *yn,
                     double *rn, double *xrn, double *yrn, double *tn,
                     double *alpn, double *alkn, double *xnew,
                     double * zxcharl, double *zycharl, double *zxchar, double *zychar,
                     double *zxchard, double *zychard, double *zxcharld,
                     double *zycharld, const std::string& typ)
{

    double tcl;
    int nchar;
    int nchard;
    int ncharl;
    int ncharld;
    
    
    
    
        int j1;
    double xx, yy, xx0, yy0, xsr, ysr;
    int irep;
    
    
    

    normtot(mw, xpodz, 0.0, 1.0);
    normtot(nw, ypodz, 0.0, 1.0);

    /*     parametry charakterystyczne linii obwiedniowej */

    double tcharld = gett(zxcharld, zycharld, *nn, xn, yn, rn, xrn, yrn, tn);
    double tcharl = gett(zxcharl, zycharl, *nn, xn, yn, rn, xrn, yrn, tn);
    double tchar = gett(zxchar, zychar, *nn, xn, yn, rn, xrn, yrn, tn);
    double tchard = gett(zxchard, zychard, *nn, xn, yn, rn, xrn, yrn, tn);

    *zxcharld = getxbrz(tcharld, tn, xn, rn, xrn, alpn, alkn);
    *zycharld = getybrz(tcharld, tn, yn, rn, yrn, alpn, alkn);
    *zxcharl = getxbrz(tcharl, tn, xn, rn, xrn, alpn, alkn);
    *zycharl = getybrz(tcharl, tn, yn, rn, yrn, alpn, alkn);
    *zxchar = getxbrz(tchar, tn, xn, rn, xrn, alpn, alkn);
    *zychar = getybrz(tchar, tn, yn, rn, yrn, alpn, alkn);
    *zxchard = getxbrz(tchard, tn, xn, rn, xrn, alpn, alkn);
    *zychard = getybrz(tchard, tn, yn, rn, yrn, alpn, alkn);
    
    double xcharld = *zxcharld;
    double ycharld = *zycharld;
    double xcharl = *zxcharl;
    double ycharl = *zycharl;
    double xchar = *zxchar;
    double ychar = *zychar;
    double xchard = *zxchard;
    double ychard = *zychard;
    

    ncharld = 1;
    while (tn[ncharld - 1] + 1e-10 < tcharld)
    {
        ++ncharld;
    }
    --ncharld;
    ncharl = 1;
    while (tn[ncharl - 1] + 1e-10 < tcharl)
    {
        ++ncharl;
    }
    --ncharl;
    nchar = 1;
    while (tn[nchar - 1] - 1e-10 < tchar)
    {
        ++nchar;
    }
    --nchar;
    nchard = 1;
    while (tn[nchard - 1] + 1e-10 < tchard)
    {
        ++nchard;
    }
    --nchard;

    for (int i = 1; i < ncharl + 1; ++i)
    {
        xn[*nn + i] = xn[i];
        yn[*nn + i] = yn[i];
        rn[*nn + i] = rn[i];
        xrn[*nn + i] = xrn[i];
        yrn[*nn + i] = yrn[i];
        alpn[*nn + i] = alpn[i];
        alkn[*nn + i] = alkn[i];
        tn[*nn + i] = tn[i];
    }
    tcl = tcharl + tn[*nn];
    *nn += ncharl;

    std::string finame;
    finame.assign(typ);
    finame.append("_m.otl");


    FM.sett(finame, GLOB.mnmx.sc, *nn, xn, yn, rn, xrn, yrn, alpn, alkn, tn);

#ifdef GPU_INCLUDED

    cl::Buffer Bypodz(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * NWMAX, ypodz);
    cl::Buffer Bxpodz(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * MWMAX, xpodz);
    cl::Buffer Bxn(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * 2 * IZARMAX, xn);
    cl::Buffer Byn(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * 2 * IZARMAX, yn);
    cl::Buffer Brn(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * 2 * IZARMAX, rn);
    cl::Buffer Bxrn(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * 2 * IZARMAX, xrn);
    cl::Buffer Byrn(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * 2 * IZARMAX, yrn);
    cl::Buffer Balpn(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * 2 * IZARMAX, alpn);
    cl::Buffer Balkn(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * 2 * IZARMAX, alkn);
    cl::Buffer Btn(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * 2 * IZARMAX, tn);
    cl::Buffer Bxnew(*GPU.context, CL_MEM_WRITE_ONLY | CL_MEM_USE_HOST_PTR, sizeof (double) * 2 * MWMAX * NWMAX, xnew);


    (*GPU.kernel_podzial).setArg(0, Bypodz);
    (*GPU.kernel_podzial).setArg(1, Bxpodz);
    (*GPU.kernel_podzial).setArg(2, Bxn);
    (*GPU.kernel_podzial).setArg(3, Byn);
    (*GPU.kernel_podzial).setArg(4, Brn);
    (*GPU.kernel_podzial).setArg(5, Bxrn);
    (*GPU.kernel_podzial).setArg(6, Byrn);
    (*GPU.kernel_podzial).setArg(7, Balpn);
    (*GPU.kernel_podzial).setArg(8, Balkn);
    (*GPU.kernel_podzial).setArg(9, Btn);
    (*GPU.kernel_podzial).setArg(10, Bxnew);
    (*GPU.kernel_podzial).setArg(11, *nn);
    (*GPU.kernel_podzial).setArg(12, tcharl);
    (*GPU.kernel_podzial).setArg(13, tchar);
    (*GPU.kernel_podzial).setArg(14, tchard);
    (*GPU.kernel_podzial).setArg(15, tcharld);
    (*GPU.kernel_podzial).setArg(16, xcharld);
    (*GPU.kernel_podzial).setArg(17, ycharld);
    (*GPU.kernel_podzial).setArg(18, xcharl);
    (*GPU.kernel_podzial).setArg(19, ycharl);
    (*GPU.kernel_podzial).setArg(20, xchard);
    (*GPU.kernel_podzial).setArg(21, ychard);
    (*GPU.kernel_podzial).setArg(22, xchar);
    (*GPU.kernel_podzial).setArg(23, ychar);
    (*GPU.kernel_podzial).setArg(24, ncharld);
    (*GPU.kernel_podzial).setArg(25, ncharl);
    (*GPU.kernel_podzial).setArg(26, nchard);
    (*GPU.kernel_podzial).setArg(27, nchar);
    (*GPU.kernel_podzial).setArg(28, tcl);
    (*GPU.kernel_podzial).setArg(29, nw);
    (*GPU.kernel_podzial).setArg(30, mw);
    
    
//cl::Event event;

    (*GPU.queue).enqueueNDRangeKernel((*GPU.kernel_podzial), cl::NullRange, cl::NDRange(mw, nw), cl::NullRange);//, NULL, &event);

    (*GPU.queue).enqueueReadBuffer(Bxnew, CL_TRUE, 0, sizeof (double) * 2 * MWMAX * NWMAX, xnew);
    
//    (*GPU.queue).finish();
//    cl_ulong end = event.getProfilingInfo<CL_PROFILING_COMMAND_END>();
//    cl_ulong start = event.getProfilingInfo<CL_PROFILING_COMMAND_START>();
//    std::cout << "OUT " << (end-start)/1000000. << std::endl;

#else

    for (int i = 0; i < mw; ++i)
    {
        for (int j = 0; j < nw; ++j) 
        {
            //j1 = i1 + j + 1;
            j1 = i * nw + j + 1;
            xx = getx_(ypodz[j], xpodz[i], *nn, xn, yn, rn,
                       xrn, yrn, xcharld, ycharld, xcharl, ycharl,
                       ncharld, ncharl, xchard, ychard, xchar, ychar, nchard, nchar);
            xx0 = 1e10;
            yy = 1e10;
            irep = 0;
            while (std::abs(xx - xx0) > 1e-12 && irep < 10)
            {
                ++irep;
                yy0 = yy;
                yy = gety(xx, ypodz[j], tcharl, tchar, tchard, tcharld, tn, xn, yn, rn, xrn, yrn, alpn, alkn);
                xx0 = xx;
                xx = getx_p(yy, xpodz[i], tchar, tchard, tcharld, tcl, tn, xn, yn, rn, xrn, yrn, alpn, alkn);
            }
            if (irep >= 10)
            {
                xsr = xx;
                irep = 0;
                while (std::abs(xx - xx0) > 1e-12 && irep < 200)
                {
                    ++irep;
                    yy0 = yy;
                    yy = gety(xsr, ypodz[j], tcharl, tchar, tchard, tcharld, tn, xn, yn, rn, xrn, yrn, alpn, alkn);
                    ysr = (yy0 + yy) / 2.;
                    xx0 = xx;
                    xx = getx_p(ysr, xpodz[i], tchar, tchard, tcharld, tcl, tn, xn, yn, rn, xrn, yrn, alpn, alkn);
                    xsr = (xx0 + xx) / 2.;
                }
            }
            xnew[2 * j1 - 2] = xx;
            xnew[2 * j1 - 1] = yy;
        }
    }
    
#endif    
    return 0;
}


#ifndef GPU_INCLUDED
double Podzial::getx_(double eta, double ksi, int nn, double *
                      xn, double *yn, double *rn, double *xrn, double *yrn,
                      double xcharld, double ycharld, double xcharl,
                      double ycharl, int ncharld, int ncharl, double
                      xchard, double ychard, double xchar, double ychar,
                      int nchard, int nchar)
{


    /*     fl=getxbel(eta,0.d0,xcharld,1.d0,xcharl,yn(2),xn(2), */
    /*    1                                           nn-ncharld+ncharl) */

    //int i__1 = nn - ncharld + ncharl;
    double fl = getxbep_(eta, ycharld, xcharld, ycharl, xcharl, &yn[ncharld - 2],
                         &xn[ncharld - 2], &rn[ncharld - 2], &yrn[ncharld - 2], &xrn[ncharld - 2]);
    //i__1 = nchard - nchar + 1;
    double fp = getxbep_(eta, ychard, xchard, ychar, xchar, &yn[nchar - 2], &xn[nchar - 2],
                         &rn[nchar - 2], &yrn[nchar - 2], &xrn[nchar - 2]);

    return fl + ksi * (fp - fl);
}

double Podzial::getx_p(double y, double ksi, double t1p,
                       double t2p, double t1l, double t2l, double *tn,
                       double *xn, double *yn, double *rn, double *xrn,
                       double *yrn, double *alpn, double *alkn)
{
    double tbiezp = getty(t1p, t2p, y, tn, yn, rn, yrn, alpn, alkn);
    double fp = getxbrz(tbiezp, tn, xn, rn, xrn, alpn, alkn);
    double tbiezl = getty(t1l, t2l, y, tn, yn, rn, yrn, alpn, alkn);
    double fl = getxbrz(tbiezl, tn, xn, rn, xrn, alpn, alkn);

    return fl + ksi * (fp - fl);
}

double Podzial::getxbep_(double eta, double xp, double yp,
                         double xk, double yk, double *xn, double *yn,
                         double *rn, double *xrn, double *yrn)
{

    double yfu;
    double wsp = xk - xp;
    double etan = xp + eta * wsp;

    if (etan >= xk)
    {
        yfu = yk;
    }
    else if (etan <= xp)
    {
        yfu = yp;
    }
    else
    {
        double dx, dy, rr, dxb, dyb, cosal, sinal;
        int i = 1;
        while (etan < xn[i])
        {
            ++i;
        }
        rr = std::abs(rn[i - 1]);
        if (rr < 1e-10)
        {
            dx = xn[i] - xn[i - 1];
            dy = yn[i] - yn[i - 1];
            dxb = etan - xn[i - 1];
            dyb = dy * dxb / dx;
            yfu = yn[i - 1] + dyb;
        }
        else
        {
            sinal = (etan - xrn[i - 1]) / rr;
            cosal = std::sqrt(1. - sinal * sinal);
            yfu = yrn[i - 1] + rr * cosal;
        }
    }

    return yfu;
}

double Podzial::gety(double x, double eta, double t1g, double
                     t2g, double t1d, double t2d, double *tn, double *xn,
                     double *yn, double *rn, double *xrn, double *yrn,
                     double *alpn, double *alkn)
{

    double tbiezg = gettx(t1g, t2g, x, tn, xn, rn, xrn, alpn, alkn);
    double fg = getybrz(tbiezg, tn, yn, rn, yrn, alpn, alkn);
    double tbiezd = gettx(t1d, t2d, x, tn, xn, rn, xrn, alpn, alkn);
    double fd = getybrz(tbiezd, tn, yn, rn, yrn, alpn, alkn);

    return fd + eta * (fg - fd);
}

double Podzial::gettx(double t1, double t2, double x, double *
                      tn, double *xn, double *rn, double *xrn, double *alpn,
                      double *alkn)
{

    double xx, tbiez;
    double tstart = t1;
    double tstop = t2;
    double xstart = getxbrz(tstart, tn, xn, rn, xrn, alpn, alkn);

    while (std::abs(tstop - tstart) > 1e-10)
    {
        tbiez = (tstart + tstop) / 2.;
        xx = getxbrz(tbiez, tn, xn, rn, xrn, alpn, alkn);
        if ((x - xx) * (xx - xstart) > 0.)
        {
            tstart = tbiez;
            xstart = xx;
        }
        else
        {
            tstop = tbiez;
        }
    }

    return (tstart + tstop) / 2.;
}

double Podzial::getty(double t1, double t2, double y, double *
                      tn, double *yn, double *rn, double *yrn, double *alpn,
                      double *alkn)
{

    double yy, tbiez;
    double tstart = t1;
    double tstop = t2;
    double ystart = getybrz(tstart, tn, yn, rn, yrn, alpn, alkn);

    while (std::abs(tstop - tstart) > 1e-10)
    {
        tbiez = (tstart + tstop) / 2.;
        yy = getybrz(tbiez, tn, yn, rn, yrn, alpn, alkn);
        if ((y - yy) * (yy - ystart) > 0.)
        {
            tstart = tbiez;
            ystart = yy;
        }
        else
        {
            tstop = tbiez;
        }
    }

    return (tstart + tstop) / 2.;
}
#endif


double Podzial::getxbrz(double tbiez, double *tn, double *xn, double *rn, double *xrn, double *alpn, double *alkn)
{
    int i = 1;

    while (tbiez > tn[i] + 1e-10)
    {
        ++i;
    }
    double tloc = tbiez - tn[i - 1];

    return xbrz(tloc, tn[i - 1], xn[i - 1], tn[i], xn[i], rn[i - 1], xrn[i - 1], alpn[i - 1], alkn[i - 1]);

}

double Podzial::getybrz(double tbiez, double *tn, double *yn,
                        double *rn, double *yrn, double *alpn, double *alkn)
{

    int i = 1;
    while (tbiez > tn[i] + 1e-10)
    {
        ++i;
    }

    double tloc = tbiez - tn[i - 1];

    return ybrz0(tloc, tn[i - 1], yn[i - 1], tn[i], yn[i], rn[i - 1], yrn[i - 1], alpn[i - 1], alkn[i - 1]);
}

double Podzial::getxbep(double eta, double xp, double yp,
                        double xk, double yk, double *xn, double *yn, int n)
{

    double xkk, ykk, xpp, ypp;

    if (eta >= xn[0])
    {
        xpp = xp;
        ypp = yp;
        xkk = xn[0];
        ykk = yn[0];
    }
    else if (eta <= xn[n])
    {
        xkk = xk;
        ykk = yk;
        xpp = xn[n];
        ypp = yn[n];
    }
    else
    {
        int i = 1;
        while (eta < xn[i])
        {
            ++i;
        }
        xpp = xn[i - 1];
        ypp = yn[i - 1];
        xkk = xn[i];
        ykk = yn[i];
    }
    double dx = xkk - xpp;
    double dy = ykk - ypp;
    double dxb = eta - xpp;
    double dyb = dy * dxb / dx;

    return ypp + dyb;
}

/* ----> poszukiwanie odcietej punktu nalezacego do brzegu sieci */

double Podzial::xbrz(double t, double t1, double x1, double t2,
                     double x2, double r, double xr, double alp, double alk)
{

    if (std::abs(r) < 1e-10)
    {
        double ax = x2 - x1;
        double tx = t2 - t1;
        return (x1 + t * ax / tx);
    }
    else
    {
        double alx = alk - alp;
        double tx = t2 - t1;
        double al = alp + t * alx / tx;
        return (xr + std::abs(r) * cos(al));
    }
}



/* ----> poszukiwanie rzednej punktu nalezacego do brzegu sieci */

/*       dla plaszczyzny wyjscia z walcow */

double Podzial::ybrz0(double t, double t1, double y1, double t2,
                      double y2, double r, double yr, double alp, double alk)
{

    if (std::abs(r) < 1e-10)
    {
        double ay = y2 - y1;
        double ty = t2 - t1;
        return (y1 + t * ay / ty);
    }
    else
    {
        double aly = alk - alp;
        double ty = t2 - t1;
        double al = alp + t * aly / ty;
        return (yr + std::abs(r) * sin(al));
    }
}

/*     szukanie w przedzialach */
double Podzial::gett(double *xp, double *yp, int n, double *x,
                     double *y, double *r, double *xr, double *yr,
                     double *t)
{
    //printf("E %f\n", *t);
    double ret_val = 0.0;
    int i;
    double v, x1, x2, y1, y2, da, db, dl, tb, ax, ay, bx, dx, dy,
            by, xl, yl, vx, vy, al1, al2, cfi, sfi, xob, yob, dxy;
    double xop, yop, ypp, xsr, ysr, xxp, yyp, dxy1, dxy2, cosfi,
            sinfi;
    int j = 0;
    double dxymin = 1e10;

    for (i = 1; i <= n; ++i)
    {
        xop = x[i - 1];
        yop = y[i - 1];
        xob = x[i];
        yob = y[i];
        if (xop > xob)
        {
            xxp = xop;
            xop = xob;
            xob = xxp;
        }
        if (yop > yob)
        {
            yyp = *yp;
            yop = yob;
            yob = yyp;
        }
        if (xop - 1e-10 < *xp && xob + 1e-10 > *xp)
        {
            xsr = (xop + xob) / 2.;
            ysr = (yop + yob) / 2.;
            dx = *xp - xsr;
            dy = *yp - ysr;
            dxy = std::sqrt(dx * dx + dy * dy);
            if (dxy < dxymin)
            {
                dxymin = dxy;
                j = i;
            }
        }
    }

    if (j > 0 && j <= n)
    {
        while (std::abs(x[i] - x[i - 1]) < 1e-10)
        {
            dx = *xp - x[i - 1];
            dy = *yp - y[i - 1];
            dxy1 = std::sqrt(dx * dx + dy * dy);
            dx = *xp - x[i];
            dy = *yp - y[i];
            dxy2 = std::sqrt(dx * dx + dy * dy);
            if (dxy1 > dxy2)
            {
                ++i;
            }
            else
            {
                --i;
            }
        }
    }
    else
    {

        /*     szukanie najblizszego przedzialu */
        double getpar;
        j = 0;
        double cfimin = 1.;

        for (int i = 1; i <= n - 1; ++i)
        {
            dx = x[i] - x[i - 1];
            dy = y[i] - y[i - 1];
            dl = std::sqrt(dx * dx + dy * dy);
            if (dl > 1e-10)
            {
                ax = x[i - 1] - *xp;
                ay = y[i - 1] - *yp;
                da = std::sqrt(ax * ax + ay * ay);
                if (da < 1e-10)
                {
                    getpar = t[i - 1];
                    return ret_val;
                }
                bx = x[i] - *xp;
                by = y[i] - *yp;
                db = std::sqrt(bx * bx + by * by);
                if (db < 1e-10)
                {
                    getpar = t[i];
                    return ret_val;
                }
                cfi = (ax * bx + ay * by) / da / db;
                if (cfi < cfimin)
                {
                    j = i;
                    cfimin = cfi;
                }
            }
        }
        x1 = x[j - 1];
        x2 = x[j];
        y1 = y[j - 1];
        y2 = y[j];
        dx = x2 - x1;
        dy = y2 - y1;
        dl = std::sqrt(dx * dx + dy * dy);
        sfi = dy / dl;
        cfi = dx / dl;
        GLOB.rot(xp, yp, sfi, cfi);
        GLOB.rot(&x1, &y1, sfi, cfi);
        GLOB.rot(&x2, &y2, sfi, cfi);
        *yp = y1;
        if (*xp < x1)
        {
            *xp = x1;
        }
        if (*xp > x2)
        {
            *xp = x2;
        }

        GLOB.rot(xp, yp, -sfi, cfi);
    }

    if (j < 1 || j > n)
    {
        FM.FL2 << std::setw(20) << *xp << std::setw(20) << *yp << "\n";

        for (int i = 0; i < n + 1; ++i)
        {
            FM.FL2 << std::setw(20) << x[i] << std::setw(20) << y[i] << std::setw(20) << r[i] << "\n";
        }
        //exit(0);
    }
    //std::cout << "przek" << std::endl;
    i = j;
    if (std::abs(r[i - 1]) < 1e-10)
    {
        x1 = x[i - 1];
        y1 = y[i - 1];
        x2 = x[i];
        y2 = y[i];
        xl = *xp;
        yl = *yp;
        vx = x2 - x1;
        vy = y2 - y1;
        v = std::sqrt(vx * vx + vy * vy);
        sinfi = vy / v;
        cosfi = vx / v;
        GLOB.rot(&x1, &y1, sinfi, cosfi);
        GLOB.rot(&x2, &y2, sinfi, cosfi);
        GLOB.rot(&xl, &yl, sinfi, cosfi);
        if (xl < x1)
        {
            xl = x1;
        }
        if (xl > x2)
        {
            xl = x2;
        }
        tb = xl - x1;

    }
    else
    {
        //std::cout << x[i - 1] << " " <<  y[i - 1] << std::endl;
        //std::cout << *xr << " " << *yr << " " << -i << std::endl;
        ypp = ypasma(*xp, n, x, y, r, xr, yr, -i);

        GLOB.getalf(x[i - 1], y[i - 1], *xp, ypp, r[i - 1], xr[i - 1], yr[i - 1], &al1, &al2, &tb);

        //std::cout << al1 << " " << al2 << std::endl;
        //        getalf(x[i], y[i], *xp, ypp, r[i], xr[i], yr[i], &al1, &al2, &tb);
    }
    //std::cout << t[i - 1] + tb << std::endl;
    return t[i - 1] + tb;
    //return t[i] + tb;
}



/* ----> podprogram ma za zadanie unormowanie tablicy podzialu (podz */
/* ----> tak, aby zachowac stosunki kolejnych odcinkow i aby ostatnia */
/* ----> wartoscia bylo tk, a pierwsza tp. */
/* ----> znaczenie parametrow: */
/*     ----> podz  - tablica podzialu */
/*     ----> n     - ilosc elementow macierzy podz */

/*     ----> tp,tk - zakres zmian tablicy podz */

int Podzial::normtot(int n, double *podz, double tp, double tk)
{

    double podzp = podz[0];
    double zakt = tk - tp;
    double zakp = podz[n - 1] - podzp;

    for (int i = 0; i < n; ++i)
    {
        podz[i] = tp + (podz[i] - podzp) * zakt / zakp;
    }
    return 0;
}


/* ----> funkcja ma za zadanie wyznaczenie wartosci wspolrzednej y-kowej */
/* ----> punktu polozonego na brzegu przekroju poprzecznego pasma wcho- */
/* ----> dzacego miedzy walce. */
/* ----> Uklad wspolrzednych jest tak usytuowany, ze poczatek ukladu */
/* ----> lezy w geometrycznym srodku wykroju walcow na plaszczyznie */
/* ----> przechodzacej przez osie walcow. Os y przebiega pionowo do */
/* ----> gory, os x w prawo, a os z prostopadle do tego przekroju i */
/* ----> przeciwnie do kierunku walcowania. */
/* ----> zmienne n,x,y,r,xr,yr maja identyczny sens jak w procedurze */
/* ----> bruzda (dla pasma). Tam tez mozna znalezc ich dokladny opis. */
/* ----> imode - wskaznik sterujacy okreslajacy zachowanie sie podpro- */
/* ---->         gramu w wypadku, gdy xp nie miesci sie w zakresie */
/* ---->         (x(1),x(n+1)): */
/*     ----> imode< 0 - punkt jest w zakresie i - i+1 */
/*     ----> imode= 0 - ypasma przyjmuje wartosc rowna: */
/*         ---->        y(1)   dla xp < x(1) */
/*         ---->        y(n+1) dla xp > x(n+1) */
/*     ----> imode= 1 - ypasma przyjmuje wartosc nieokreslana i nastepuje */

/*     ---->     wstrzymanie obliczen z odpowiednim komunikatem */


double Podzial::ypasma(double xp, int n, double *x, double *y,
                       double *r, double *xr, double *yr, int imode)
{
    //std::cout << "YPAS" << std::endl;
    double ret_val;
    int i = 0;

    if (imode < 0)
    {
        i = -imode;
    }
    else
    {
        bool go_l2 = false;
        for (i = 1; i <= n; ++i)
        {
            if (xp >= x[i - 1] && xp <= x[i])
            {
                go_l2 = true;
                break;
            }
        }

        if (!go_l2)
        {
            if (imode == 0)
            {
                if (xp < x[0])
                {
                    ret_val = y[0];
                }
                else
                {
                    ret_val = y[n];
                }
                return ret_val;
            }
            else
            {
                std::cout << " Bledy przy wyznaczaniu punktu na powierzchni pasma" <<
                        "\nWspolrzedna x-owa poczatku wykroju    :  " << *x <<
                        "\nWspolrzedna x-owa konca wykroju      :  " << x[n] <<
                        "\nZadana wspolrzedna x-owa punktu      :  " << xp << "\n";
            }
        }
    }

    double x1 = x[i - 1];
    double y1 = y[i - 1];
    double x2 = x[i];
    double y2 = y[i];

    if (std::abs(r[i - 1]) > 0.0)
    {
        double ro = r[i - 1];
        double xro = xr[i - 1];
        double yro = yr[i - 1];
        double rx = xp - xro;
        double rx2 = ro * ro - rx * rx;
        if (rx2 >= 0.)
        {
            ret_val = yro + ro * std::sqrt(rx2) / std::abs(ro);
        }
        else
        {
            ret_val = yro;
        }
    }
    else
    {
        ret_val = y1 + (y2 - y1) * (xp - x1) / (x2 - x1);
    }

    return ret_val;

}
